Ce package implémente la méthode du graphical lasso avec pénalité l1. Il permet d’estimer une matrice de précision creuse, qui peut correspondre à une variable suivant une normale multi-variée. La matrice de précision donne dans ce cas accès aux arrêtes du graph représentant les relations d’indépendances conditionnelles au sein de la variable.
Ci-dessous on explicite l’exemple donné dans le package.
set.seed(100)
x<-matrix(rnorm(50*20),ncol=20)
s<- var(x)
On considère une matrice \(X\) de taille 50*20, correspondant à 50 échantillons d’une variable de dimension 20. Cette variable suit une loi normale multi-variée centrée, de variance \(s\).
glasso prend en entrée à minima la matrice de covariance et le paramètre de pénalisation rho (rho=0 signifie pas de régularisation). Il est également possible de forcer certains coefficients de l’estimée de l’inverse de la matrice de covariance à zéro, afin d’estimer un graph avec des arêtes manquantes.
En sortie on récupère une liste de 7 éléments
a<-glasso(s, rho=.01)
aa<-glasso(s,rho=.1, w.init=a$w, wi.init=a$wi)
Un autre exemple : introduction de zéros structurels et absence de régularisation (Whittaker’s book on graphical models) :
S=matrix(0,nrow=4,ncol=4)
s=c(10,1,5,4,10,2,6,10,3,10)
S[row(S)>=col(S)]=s # remplissage de la partie triangulaire inférieure
S=(S+t(S))
diag(S)<-10 # symmétrisation en gardant la diagonale
zero<-matrix(c(1,3,2,4),ncol=2,byrow=TRUE)
a_zeros<-glasso(S,0,zero=zero)
## Warning in glasso(S, 0, zero = zero): With rho=0, there may be convergence problems if the input matrix
## is not of full rank
Dans cet exemple, on a défini la matrice des zéros. Lignes en colonne 1, colonne en colonne deux. On observe bien des 0 aux positions [1,3] et [2,4] de la matrice wi.
library(ggplot2)
library(reshape2)
library(grid)
library(gridExtra)
library(RColorBrewer)
pal<-brewer.pal(8, "Spectral")
g1<-ggplot(data = melt(a$w), aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+
labs(title="K, rho = 0.01")+
scale_fill_gradient2(low="red", high="blue",mid="white",midpoint=0)
g2<-ggplot(data = melt(a$wi), aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+
labs(title="inverse de K, rho = 0.01")+
scale_fill_gradient2(low="red", high="blue",mid="white",midpoint=0)
g3<-ggplot(data = melt(aa$w), aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+
labs(title="K à partir d'une première estimation,\n rho = 0.1")+
scale_fill_gradient2(low="red", high="blue",mid="white",midpoint=0)
g4<-ggplot(data = melt(aa$wi), aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+
labs(title="inverse de K à partir d'une première estimation,\n rho = 0.1")+
scale_fill_gradient2(low="red", high="blue",mid="white",midpoint=0)
g5<-ggplot(data = melt(a_zeros$wi), aes(x=Var1, y=Var2, fill=value)) +
geom_tile()+
labs(title="inverse, zéros forcés, rho=0")+
scale_fill_gradient2(low="red", high="blue",mid="white",midpoint=0)
grid.arrange(g1, g2, g3, g4,g5)
Tout ce qui suit vient du site génial de Katherine Ognyanova <www.kateto.net/network-visualization>. Les packages de base sont igraph, network, sna, visNetwork, threejs,, networkD3 et ndtv. A noter le bug du package igraph, résolu sur github.
library(network)
library(visNetwork)
library(sna)
library(ndtv)
# install.packages("devtools")
# library(devtools)
# install_github("igraph/rigraph")
library(igraph) #bug on the CRAN package, resolved on github
library(threejs)
library(networkD3)
Pour la suite nous allons continuer avec les deux exemples, afin d’avoir différentes dimensions de réseaux.
La matrice d’adjacence est composée de 0 et de valeurs non nulles qui correspondent aux indicatrices de liens entre les différents noeuds du réseau. Dans nos exemple on reste dans le cas gaussien : on va seuiller la matrice inverse de covariance, estimée par glasso, afin de déduire les arêtes du graph.
Le package igraph permet de créer des graphs à partir de différents format : data frame, matrice d’adjacence ou matrice d’incidence. Par data frame, la fonction va prendre en paramètres deux sources de données : les liens et les noeuds, dont chacun peut avoir plusieurs autres colonnes renseignant des attributs supplémentaires. Le package peut aussi tout construire à partir de la matrice d’adjacence (matrice carrée, pas forcément symmétrique si dirigé, et les valeurs peuvent signifier des poids si elles sont différentes de 1). La matrice d’incidence enfin n’est pas forcément carrée et représente les liens entre deux sets de variables, et les sets ne peuvent pas avoir de liaisons internes.
On a les fonctions :
library(tidyr)
library(dplyr)
library(magrittr)
matrice_adj<-function(precision,seuil){
adj<-unlist(apply(precision,1, function(x) ifelse(abs(x)<seuil,0,1)))
return(adj)
}
build_net<-function(precision,seuil, multiple, boucles){
matrice<-data.frame(matrice_adj(precision,seuil))
matrice$start<-colnames(matrice)
links<-filter(gather(matrice,vertices,status,-start),status==1)[,1:2]
nodes<-colnames(matrice)[-grep("start",colnames(matrice))]
net <- graph_from_data_frame(d=links, vertices=nodes, directed=F)%>% simplify( remove.multiple = multiple, remove.loops = boucles)
return(net)
}
net_from_matrix<-function(precision,seuil, boucles){
matrice<-matrice_adj(precision,seuil)
net <- graph_from_adjacency_matrix(matrice, mode="upper", diag = boucles)
return(net)
}
net_adj<-net_from_matrix(aa$wi,0.1,F)
net2_adj<-net_from_matrix(a_zeros$wi,0.005,F)
net<-build_net(aa$wi,0.1,T,T)
net2<-build_net(a_zeros$wi,0.005,T,F)
l1 <- layout_with_fr(net)
l1 <- norm_coords(l1, ymin=-1, ymax=1, xmin=-1, xmax=1)
l2 <- layout_with_fr(net_adj)
l2 <- norm_coords(l2, ymin=-1, ymax=1, xmin=-1, xmax=1)
l3 <- layout_with_fr(net2)
l3 <- norm_coords(l3, ymin=-1, ymax=1, xmin=-1, xmax=1)
l4 <- layout_with_fr(net2_adj)
l4 <- norm_coords(l4, ymin=-1, ymax=1, xmin=-1, xmax=1)
par(mfrow=c(2,2))
plot(net,main="dim=20", rescale=F, layout=l1,vertex.label=NA,margin=c(0,0,0,0))
plot(net_adj,main="dim 20 adj", rescale=F, layout=l2,margin=c(0,0,0,0))
plot(net2, main="Whittaker", rescale=F, layout=l3,margin=c(0,0,0,0))
plot(net2_adj, main="Whittaker adj", rescale=F, layout=l4,margin=c(0,0,0,0))
deg <- degree(net_adj, mode="out")
V(net_adj)$size <- deg*2+5
V(net_adj)$color=pal
## Warning in vattrs[[name]][index] <- value: le nombre d'objets à remplacer
## n'est pas multiple de la taille du remplacement
E(net_adj)$color=pal[7]
V(net_adj)$label = NA
E(net_adj)$curved=.1
V(net_adj)$frame.color=pal[1]
plot(net_adj,margin=c(0,0,0,0) )
plot(net_adj, layout=layout_randomly)
Calcul des coordonnées du réseau à l’avance :
l <- layout_in_circle(net_adj)
plot(net_adj, layout=l)
l <- layout_on_sphere(net_adj)
plot(net_adj, layout=l)
L’algorithme de Fruchterman-Reingold est très utilisé. Il optimise la disposition des points : arêtes de taille similaires et minimum de croisements. Le résultat n’est pas fixe sur plusieurs simulations. Pour garder la même disposition, on peut sauver le layout dans une variable l.
L’algorithme de Kamada Kawai est similaire : il tente d’optimier la structure du réseau dans l’espace. L’algorithme LGL permet de spécifier une racine au centre du graph. Enfin l’algorithme MDS permet de spécifier une matrice de distance, ce qui donne une interprétation claire des positions et des distances entre noeuds. Cependant les noeuds peuvent se superposer ce qui limite la clarté du graph. L’algorithme essaie de repérer des similarités entre les noeuds.
l <- layout_with_fr(net_adj)
plot(net_adj, layout=l)
l <- layout_with_kk(net_adj)
plot(net_adj, layout=l)
plot(net_adj, layout=layout_with_lgl, root=5)
plot(net_adj, layout=layout_with_mds)
cut.off <- mean(deg)
net.sp <- delete_edges(net_adj, E(net_adj)[deg<cut.off])
plot(net.sp) #je ne comprends pas ce qu'il fait, ce code est plus pour la démo
Le package igraph implémente une large palette de méthodes pour la clusterisation du grpah :
# Community detection returns an object of class "communities"
# which igraph knows how to plot:
clust<-function(methode){
if(methode=="cluster_spinglass"){
net<-delete_vertices(net_adj,deg<1)
}else{
net<-net_adj
}
clp <- get(methode)(net)
plot(clp, net, main=methode)
}
par(mfrow=c(5,2))
clust("cluster_edge_betweenness")
clust("cluster_fast_greedy")
clust("cluster_infomap")
clust("cluster_label_prop")
clust("cluster_leading_eigen")
clust("cluster_louvain")
clust("cluster_optimal")
clust("cluster_spinglass")
clust("cluster_walktrap")
# We can also plot the communities without relying on their built-in plot:
clp <- cluster_walktrap(net)
V(net_adj)$community <- clp$membership
plot(net_adj, vertex.color=pal[V(net_adj)$community])
Ce code permet de colorier le graph selon la longueur de propagation depuis le noeud numéro 3. Il faut noter que ce code ne fonctionne pas en présence de noeud isolé (Inf dans les distances).
net<-delete_vertices(net_adj,deg<1)
dist <- distances(net, v=V(net)[3],
to=V(net), weights=NA)
# Set colors to plot the distances:
oranges <- colorRampPalette(c("dark red", "gold"))
col <- oranges(max(dist))
col <- col[dist+1]
plot(net,vertex.color=col, vertex.label = dist)
On peut également mettre en valeur un chemin particulier. '{A} noter les fonctions shortest_paths, ecount, ew, vcount et les attributs vpath et epath.
news.path <- shortest_paths(net_adj,
from = V(net_adj)[5],
to = V(net_adj)[6],
output = "both") # both path nodes and edges
# Generate edge color variable to plot the path:
ecol <- rep("gray80", ecount(net_adj))
ecol[unlist(news.path$epath)] <- "orange"
# Generate edge width variable to plot the path:
ew <- rep(2, ecount(net_adj))
ew[unlist(news.path$epath)] <- 4
# Generate node color variable to plot the path:
vcol <- rep("gray40", vcount(net_adj))
vcol[unlist(news.path$vpath)] <- "gold"
plot(net_adj, vertex.color=vcol, edge.color=ecol,
edge.width=ew, edge.arrow.mode=0 )
Ici on utilise les fonction incident et neighbors.
nb<-12
inc.edges <- incident(net_adj, V(net_adj)[nb], mode="all")
# Set colors to plot the selected edges.
ecol <- rep("gray80", ecount(net_adj))
ecol[inc.edges] <- "orange"
ew <- rep(2, ecount(net_adj))
ew[inc.edges] <- 4
vcol <- rep("grey40", vcount(net_adj))
vcol[V(net_adj)[nb]] <- "gold"
neigh.nodes <- neighbors(net_adj, V(net_adj)[12], mode="out")
# Set colors to plot the neighbors:
vcol[neigh.nodes] <- "#ff9d00"
plot(net_adj, vertex.color=vcol, edge.color=ecol,edge.width=ew)
Il est aussi possible de “marquer” des groupes.
plot(net_adj, mark.groups=list(c(1 , 9, 10 ,15 ,16 ,19), c(4,6,7)),
mark.col=c("#C5E5E7","#ECD89A"), mark.border=NA)
Pour déplacer les noeuds à l’envie. On récupère les coordonnées modifiées à la main.
tkid <- tkplot(net_adj) #tkid is the id of the tkplot that will open
l <- tkplot.getcoords(tkid) # grab the coordinates from tkplot
plot(net_adj, layout=l)
## [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
## [18] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
Un réseau en utiilisant le package visNetwork. Besoin des dataframes de noeuds et d’arêtes.
v<-visIgraph(net_adj,physics=TRUE)
visClusteringByHubsize(v)
NetworkD3
matrice<-data.frame(matrice_adj(aa$wi,0.1))
matrice$start<-colnames(matrice)
links<-filter(gather(matrice,vertices,status,-start),status==1)[,1:2]
colnames(links)<-c("from","to")
el <- data.frame(from=as.numeric(factor(links$from))-1,
to=as.numeric(factor(links$to))-1 )
nodes<-data.frame(el$from[which(!duplicated(el$from))])
colnames(nodes)<-c("id")
nodes$level<-rep(seq(1:4),each=5)
nl <- cbind(idn=factor(nodes$id, levels=nodes$id), nodes)
forceNetwork(Links = el, Nodes = nl, Source="from", Target="to",
NodeID = "idn", Group = "level",linkWidth = 1,
linkColour = "#afafaf", fontSize=12, zoom=F, legend=T,
Nodesize=2, opacity = 0.8, charge=-20,
width = 600, height = 400)